Function RICHARDSON_LUCY

Package

Short description

Usage

    y = richardson_lucy(z, filter, options)

Input

    z (numeric): The filtered signal to be deconvolved.
    filter (struct): The filter used to obtain z, in the format of filters
        from FILTER_BANK.
    options (struct, optional): Different parameters, such as
            rl_iter (numeric): The iterations to output. If multiple are
                specified, they will be arranged in y as columns (default 32)

Output

     y (numeric): The result of the algorithm.

Description

    The Richardson-Lucy algorithm estimates a signal y from its convolution
    z = y*h, where h is a filter, often lowpass. The algorithm is iterative
    and converges to the pseudo-inverse of the convolution operator, which is
    unstable when the Fourier transform of the filter takes values close to
    zero. As a result, only a few iterations are used. For more details, see
    [1] and [2].
 References
    [1] L. Lucy, “An iterative technique for the rectification of observed 
       distributions,” Astron. J., vol. 79, p. 745, 1974. 
    [2] W.H. Richardson, "Bayesian-based iterative method of image 
       restoration." JOSA vol. 62.1, p. 55-59, 1972.
function yn = richardson_lucy(z, filter, options)
	if nargin < 3
		options = struct();
	options = fill_struct(options, 'rl_iter', 32);
	yn = zeros(length(z),length(options.rl_iter));
	h = realize_filter(filter, size(z,1));
	h_filter = @(u)(real(ifft(fft(u).*h)));
	ht_filter = @(u)(real(ifft(fft(u).*conj(h))));
	y = z;
	for n = 1:max(options.rl_iter)
		zr = z./(h_filter(y));
		y = y.*ht_filter(zr);
		if ismember(n, options.rl_iter)
			[~,k] = ismember(n, options.rl_iter);
			yn(:,k) = y;

See also